NUMERICAL STUDY OF THE DERIVATIVE OF THE 
RIEMANN ZETA FUNCTION AT ZEROS 
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Dedicated to Professor Akio Fujii on his retirement. 

Abstract. The derivative of the Riemann zeta function was computed nu- 
merically on several large sets of zeros at large heights. Comparisons to known 
and conjectured asymptotics are presented. 

1. Introduction 

Throughout this paper, we assume the truth of the Riemann Hypothesis (RH), 
and we let 7 „ > denote the ordinate of the n-th non-trivial zero of ((s). Hejhal [Hej 
assumed the RH and a weak consequence of Montgomery's |Moj pair-correlation 
conjecture, namely that for some r > 0, there is a constant B such that 

(1.1) limsup — |{n : N<n<2N, ( 7n+1 - 7n ) log 7n < c}| < Bc T , 

N->oo J v 

holds for all c € (0, 1). Under these assumptions, he proved the following central 

limit theorem: for a < j3, 

(1.2) 

27rC'(l/2 + i7 „) 



lim 



1 

N 



log 



N < n < 2N . 



log( 7 „/27r) 



\ log log N 



€ {a J) 



1 



e- x2 / 2 dx. 



So under these assumptions, log \('(l/2 + « 7n )|, suitably normalized, converges in 
distribution over fixed ranges to a standard normal variable. To obtain more precise 
information about the tails of the distribution, we consider the moments 



(1.3) 



Jx(T) 



— y 

N(T) ^ 

V ; 0<7„<T 



|C'(l/2- 



\2\ 



where N(T) :— J2o<-y n <T 1 — log ^5 + O(logT) is the zero counting function. 
Notice that J\(T) is defined for all A provided the zeros of £(s) are simple, as is 
widely believed. 

Gonek |Golj |Go2j carried out an extensive study of J\(T). He proved, under 
the assumption of the RH, that Ji(T) ~ ^(logT) 3 as T —> 00. It was suggested 
by Gonek |Go2l . and independently by Hejhal [Hej . that J\(T) is on the order of 
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(logT) A ( A+2 ). Ng [Ng proved, under the RH, that J 2 (T) is order of (logT) 8 , which 
is in agreement with that suggestion. 

Hughes, Keating, and O'Connell [HKO] . applied the random matrix philosophy 
(e.g. see |KSj ). which predicts that certain behaviors of L- functions are mimicked 
statistically by characteristic polynomials of large matrices from the classical com- 
pact groups. This led them to predict that for Re(A) > —3/2, 

(1.4) ^m-WG^O ^) as T oo, 

where G(z) is the Barnes G-function, and a(k) is an "arithmetic factor." The 



conjecture (1.4 1 is consistent with previous theorems and conjectures. 

Recently, Conrey and Snaith |CSj . assuming the ratios conjecture, gave lower 
order terms in asymptotic expansions for J\(T) and J2{T). They conjectured the 
existence of certain polynomials P\(x), for 2A = 2 and 2A = 4, such that 

(1-5) £ |C'(l/2 + * 7 „)| 2A ~ f T p x (logf) dt, 

0<7„<T ^° ^ 71 ' 

The conjecture for the case 2 A = 2 was subsequently proved by Milinovich [Mij . 
assuming the RH. It is expected that such polynomials exist for other integer values 
of A > as well. 

The purpose of this article is to study numerically various statistics of the de- 
rivative of the zeta function at its zeros. In particular, we consider the distribution 
of log |£'(l/2 + i7 n )|, moments of |£'(l/2 + i7„)|, and correlations among moments. 
The goal is to obtain more detailed information about the derivative at zeros, and 
to enable comparison with various conjectured and known asymptotics. Our com- 
putations rely on large sets of zeros at large heights that are described in detail in 

We find that the empirical distribution of log |£'(l/2 + i7„)|, normalized to have 
mean zero and standard deviation one, agrees generally well with the limiting nor- 
mal distribution proved by Hejhal, as shown in Figure [l] But the empirical mean 
and standard deviation pre-normalization are noticeably different from predicted 
ones. Also, as shown in Figure [2] the frequency of very small normalized values of 
log |C'(i/2 + «7 n )| is higher than predicted by a standard normal distribution, while 
the frequency of very large normalized values is lower than predicted. Since these 
differences appear to decrease steadily with height, however, they are probably not 
significant. 

To examine the tails of the distribution of log |£'(l/2-H 7 „)|, we present data for 
the moments of \('(l/2 + i-f n )\ over short ranges: 

(1-6) Jx(T,g):= 1 £ |C'(l/2 + ^r. 

v 1 v ' T<j n <T+H 

For large A, the empirical values of J\(T,H) deviate substantially from the values 



suggested by the leading term prediction (1.4). This is not surprising. Because for 



A large relative to T, the contribution of lower order terms is likely to dominate, 
and so the leading term asymptotic on its own may not suffice. Furthermore, the 
said deviations decrease steadily with height and they occur in a generally uniform 
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way for roughly 2A < 6, so they are consistent with the effect of "lower order terms" 
still being felt even at such relatively large heights. 

In the specific cases of the second and fourth moments of |£'(l/2 + i"f n )\: the 
conjectures of Conrey and Snaith |CS| supply lower order terms, and the agreement 
with the data is much better, as shown in Table [4] Q 

As A increases, the observed variability in the moments of |£'(l/2 + ij n )\ is 
more extreme, but it is still significantly less than we previously encountered in 
the moments of \((l/2 + it)\ (see [HO]). To illustrate, our computations of the 
twelfth moment of |£'(l/2 + ij n )\ over 15 separate sets of ps 10 9 zeros each (near 
the 10 23 -rd zero) show that the ratio of highest to lowest moment among the 15 
twelfth moments thus obtained was 2.36. In contrast, that ratio for the twelfth 
moment of |C(l/2 + it)\ was 16.34, which is significantly larger (see |HOJ). 

In general, the variability in statistical data for |£'(l/2 + ry n )| is considerably 
less than the variability in statistical data for |C(l/2 + it)\. It is not immediately 
clear why this should be so, considering, for instance, that the central limit theorem 
for log |C'(l/2 + ij n )\ is only conditional, while that for log |C(l/2 + it)\ is not, and 
both theorems scale by the same asymptotic variance. 

In the case of negative moments, our data is in agreement with Gonek's con- 
jecture ( [GoT] ) J-i(T) ~ ^(logT/(2 7 r))- 1 as T -> 00. But starting at 2A = -3, 
and as A decreases, the empirical behavior of negative moments becomes rapidly 
more erratic. For example, using the same 15 zero sets near the 10 23 -rd zero 
mentioned previously, the ratio of highest to lowest negative moment among them 
gets very large as A decreases; we obtain: 1.03, 8.45, 178.49, and 17240.99, for 
2A = —2, —3, —4, and —6, respectively (this can be deduced easily from Table [6]). 
Notice that the point 2A = —3 is special because it is where the leading term 



prediction (1.4 1 first breaks down due to a pole of order 1 in the ratio of Barnes 
G-functions. 

Extreme values of negative moments are caused by very few zeros. When 2A = 
—3, for instance, the largest observed moment among our 15 sets is 0.178047. About 
87% of this value is contributed by 4 zeros where |£'(l/2 + ?7„)| is small and equal 
to 0.002439, 0.002453, 0.004388, and 0.0043650 Such small values of |C'(1/ 2 + 
typically occur at pairs of consecutive zeros that are close to each other. For 
example, the values 0.002439 and 0.002453 occur at the following two consecutive 
zero ordinates: 



1.30664344087942265202071895041619 x 10 22 , 

( L7 ) 

1.30664344087942265202071898265199 x 10 22 . 

The above pair of zeros is separated by 0.00032, which is about 1/400 times the 
average spacing of zeros at that height (which is s» 0.128). 



It might be worth mentioning that we attempted to calculate the coefficients of lower order 
terms in the |CSj conjectures by calculating J\(T) for sufficiently many values of T, then solving 
the resulting system of equations. However, this did not yield good approximations of the coef- 
ficients (even for small A), which is not surprising, since the scale is logarithmic and the Conrey 
and Snaith expansion is only asymptotic. 

2 We checked such small values of [£'(1/2 + i-y) | by computing them in two ways, using the 
Odlyzko-Schonhagc algorithm, and using the straightforward Riemann-Siegel formula; the results 
from the two methods agreed to within ±10 — 6 
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To investigate possible correlations among values of |C'(l/2 + 27„)| 2A , we studied 
numerically the (shifted moment) function: 



(1.8) 



S\(T,H,m) :-- 



E 

T<7„<T+ff 



|C'(l/2 + n n )C , (l/2 + i7„+ m )| 



2 A 



We plotted S\(T, H,m), for several choices of A, T, and H, and as m varies. The 
resulting plots indicate there are long-range correlations among the values of the 
derivative at zeros. Unexpectedly, the tail of S-ziT,H,m) (Figure [3j right plot) 
strongly resembles the tail for the shifted fourth moment of |£(l/2 + it)| (Figure 4 
in [HO] ). 

To better understand these correlations, we considered the "spectrum" of log |£'(l/2+ 
«7n)|; see ( |2.6[ ) for a definition. A plot of the spectrum reveals sharp spikes, shown in 
Figure[5] These spikes can be explained heuristically by applying techniques already 
used by Fujii |Fu ( IFu2j and Gonek |Golj to estimate sums involving £'(1/2 + i"f n ). 

2. Numerical results 



Conjecture (1.2 1 suggests the mean and standard deviation of log|£'(l/2 + 
for zeros from near T — 1.3066434 x 10 22 (i.e. near the 10 23 -rd zero) are about 2.0 
and 1.4, respectively. This is far from the empirical mean and standard deviations 
listed in Table [l] which are 3.4907 and 1.0977. f\ Since these quantities grow very 
slowly (like log log T), these differences are probably not significant. 

TABLE 1. Summary statistics for log [£'(1/2 + *7n)| using sets of 10 7 zeros 
from different heights The column "Zero" lists the zero number near which the 
set is located. SD stands for standard deviation. 



Zero 


Min 


Max 


Mean 


SD 


10 le 


-3.7371 


7.3920 


3.1211 


1.0135 


10 20 


-3.2181 


8.0085 


3.3458 


1.0653 


10 23 


-2.9602 


8.2836 


3.4907 


1.0977 



We normalize the sequence {log|£'(l/2 + ij n )\ : N < n < N + 10 7 }, where 
N « 10 23 , to have mean zero and variance one. The distribution of the normalized 
sequence is illustrated in Figure [T] which contains two plots, one of the empirical 
density function, and another of the difference between the empirical density and 
the predicted (standard Gaussian) density ^k^e~ x ^ 2 ■ The fit in the first plot is 
visibly good, but there is a slight shift to the right about the center. This shift 
is made more visible in the second plot, which shows that the empirical density is 
generally larger than expected for x > 0, and is smaller than expected for x < 0. 

Near the tails, however, the situation is reversed. Figure [2] shows there is a 
deficiency in the occurrence of very large values of \('(l/2 + ij n )\, and an abundance 
in the occurrence of very small values. For instance, conjecture ( |1.2[ ) suggests that 
about 0.1462% of the values of \('(l/2 + ij n )\ near the 10 23 -rd zero should satisfy 



3 The mean and standard deviations listed in Table ^ change very little across different zero 
sets near the same height. For example, using a different set of 10 s zeros near the 10 23 -rd zero, 
the empirical mean is 3.4907 and the empirical standard deviation is 1.0978, which are very close 
the numbers listed in Table^ We note that the empirical mean and standard deviation are closer 
to the values suggested by the central limit theorem for characteristic polynomials of unitary 
matrices (see |HKOh , which are 3.47 and 1.12. 
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FIGURE 1. Empirical density of log |£'(l/2 + i"tn)\, after being normalized 
to have mean and standard deviation 1, using 10 7 values of log |f'(l/2 + i"f n )\ 
from near the 10 23 -rd zero (the bin size is 0.0512). The density of a standard 
normal variable (continuous line) is drawn to facilitate comparison. The right 
plot shows the diference, the empirical density minus the normal. 



|£'(l/2 + ij n )\ > 860, which is noticeably larger than the observed 0.1056%. The 
conjecture also suggests about 0.0736% of the values should satisfy \('(l/2 + ij n )\ < 
1, which is smaller than the observed 0.1051%. 

We remark the behavior near the tails becomes more consistent with expectation 
as height increase . For example, only 0.0025% of the time do we have log |£'(l/2 + 
ijn)\ > 3.2 near the 10 16 -th zero, which is far from the expected 0.068%, but the 
percentage increases to 0.040% near the 10 23 -rd zero. 



FIGURE 2. Distribution at the tails using 1.5 X 10 10 zeros near the 10 23 -rd 
zero (bin size is 0.01). 






For another measure of the quality of the fit to the standard Gaussian in Figure[T] 
we compare moments of both distributions. Table [2] shows the first few moments 
(the even moments in particular) agree reasonably well. Notice the odd moments 
tend to be negative, which is likely due to the aforementioned bias in the frequency 
of very small and very large values. 

To better understand the tails of the distribution of log \ (' (l/2+ij n )\, we consider 
the moments J\{T) defined in (1.3 1. Since we are interested in the asymptotic 



6 



G.A. HIARY AND A.M. ODLYZKO 



TABLE 2. Moments of log | £'(1/2 + i"f n )\, after being normalized to have 
mean zero and variance one, calculated using 10 zeros from near the 10 23 -rd 
zero. The third column is the moment of a standard Gaussian. 



Moment 


Derivatives 


Gaussian 


3rd 


-0.02728 





4th 


3.01364 


3 


5th 


-0.49120 





6th 


15.3053 


15 


7th 


-7.43073 





8th 


112.013 


105 


9th 


-118.588 





10th 


1116.64 


945 



behavior of J\(T), we compare against the leading term prediction (1.4) 
calculated ratios of the form 



We 



(2.1) 



l 

\B\ 



E 7eB IC'(l/2 + z 7 )l 



2 A 



a(A) 



G 2 (A+2) 



log 



A(A+2) 



G(2A+3) 

where B is a block of consecutive zeros, \B\ denotes the number of zeros in B, and 
T is the height where block B lies. If T is large enough, one expects the value of 
(2.1 ) to approach 1 as the block size \B\ increases. Table [3j which uses blocks of size 
\B\ R3 10 9 (except for the first set, which uses the first 10 s zeta zeros), shows that 
the empirical moments are significantly larger than the corresponding predictions, 
even for low moments. For example, the empirical second moments (2A = 2) near 
the 10 23 -rd zero are generally off from expectation by about 9.6%. 



Nevertheless, the ratios (2.1) appear to decrease towards the expected 1 as the 
height increases, and there is relatively little variation in the moment data for sets 
from near the same height when 2 A < 6. Both of these observations are consistent 
with the "lower order terms" still contributing significantly. 

The full moment prediction of [CS , which takes lower order terms into account, 
might lead one to expect that for 2A = 2, 2A = 4, as T — > oo, and for blocks B not 
too small compared to T, 



(2.2) 



E 

~f£B 



lC(V2 + i7)l 



2A 



PxQog(t/2ir))dt : 



where P\(x) is as given in [CS] , and f B is short for integrating over the interval 
spanned by the block B. To test this, we calculated ratios of the form 



£ 7eB IC'(l/2 + *7)l 



2A 



f B Px(log(t/2n))dt 



(2.3) 



As the block size increases, we expect ( |2.3[ ) to be significantly closer to 1 than ( |2.l[ ) 
since it relies on a more accurate prediction. This is indeed what Table[4]illustrates, 
where we see the fit to moment data is much better than we found in Table 
(We point out that in the case 2A = 4 only the first three terms in the full moment 



Noti ce if T is large compared to the length of the interval spanned by block B, the denominator 
in ratio ( 2.3 I is largely a function of T multiplied by the length of the interval spanned by B. 
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TABLE 3. Ratio j2.l| calculated with \B\ ~ 10 9 , except for the first set, 
which uses the first 10 s zeros. The column "Zero" lists the approximate zero 
number near which block B is located. 



Zero 


2A = 2 


2A = 4 


2A = 6 


2A = 8 


2A = 10 


2A = 12 


10° 


1 


1247 


3.1579 


91.856 


78341 


4.1016 


X 


10 9 


2.3478 


X 


10 16 


10 D 


i 


1 A O/l 


2.2087 


17.686 


1266.9 


l.ouo t 


X 


i n6 
1U 


A OKOQ 

4.yozo 


X 


into 


10 20 


1 


1123 


1.9102 


10.943 


422.72 


1.9904 


X 


10 5 


1.8362 


X 


10 9 


10 23 


1 


0964 


1.7645 


8.4406 


233.63 


6.4583 


X 


10 4 


2.7127 


X 


10 8 




1 


0964 


1.7603 


8.1602 


199.18 


4.1647 


X 


10 4 


1.1369 


X 


10 8 




1 


0964 


1.7598 


8.1879 


202.40 


4.3355 


X 


10 4 


1.2325 


X 


10 8 




1 


0964 


1.7629 


8.3221 


217.58 


5.2539 


X 


10 4 


1.7809 


X 


10 8 




1 


0964 


1.7630 


8.3861 


228.51 


6.2549 


X 


10 4 


2.6614 


X 


10 8 




1 


0964 


1.7600 


8.2022 


206.36 


4.6423 


X 


10 4 


1.4200 


X 


10 8 




1 


0965 


1.7642 


8.3321 


218.38 


5.3663 


X 


10 4 


1.8923 


X 


10 8 




1 


0965 


1.7612 


8.1862 


201.43 


4.3256 


X 


10 4 


1.2547 


X 


10 8 




1 


0963 


1.7590 


8.2176 


209.97 


4.8853 


X 


10 4 


1.5596 


X 


10 8 




1 


0964 


1.7654 


8.3856 


217.09 


4.8781 


X 


10 4 


1.4148 


X 


10 8 




1 


0963 


1.7616 


8.3009 


218.92 


5.4691 


X 


10 4 


1.9491 


X 


10 8 




1 


0964 


1.7585 


8.1576 


204.55 


4.6872 


X 


10 4 


1.5134 


X 


10 8 




1 


0965 


1.7615 


8.2380 


209.26 


4.7946 


X 


10 4 


1.5078 


X 


10 8 




1 


0963 


1.7586 


8.1764 


203.00 


4.4241 


X 


10 4 


1.2904 


X 


10 8 




1 


0964 


1.7603 


8.2037 


208.39 


4.9019 


X 


10 4 


1.6822 


X 


10 s 



conjecture were used, because these were the only terms provided explicitly in |CS| . 
It is likely the fit to the data will be even better if the missing terms are included.) 

TABLE 4. Ratio j2.3| calculated with \B\ S3 10 9 , except for the first set, 
which uses the first 10 s zeros. The column "Zero" lists the approximate zero 
number near which block B is located. 



Zero 


2A = 2 


2A = 4 


10 8 


1.0000 


1.0924 


10 16 


1.0000 


1.0144 


10 20 


1.0000 


1.0087 


10 23 


1.0000 


1.0074 




1.0000 


1.0050 




0.9999 


1.0047 




1.0000 


1.0064 




0.9999 


1.0065 




0.9999 


1.0048 




1.0000 


1.0072 




1.0000 


1.0055 




0.9998 


1.0042 




1.0000 


1.0079 




0.9999 


1.0057 




0.9999 


1.0039 




1.0000 


1.0057 




0.9999 


1.0040 




1.0000 


1.0049 



We remark that the five largest values of |£'(l/2 + i^ n )\ in our data set are 
« 7057, 6907, 6658, 6636, and 6399. The cumulative contribution of these large 
values to the 2A-th moment, as a percentage of the overall 2A-th moment, is listed 
in Table [5] for several A. 
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TABLE 5. Cumulative contribution percentage of the 5 largest values of 
|£'(l/2 + «7n)| to the empirical 2A-th moment for 1.5 X 10 10 zeros near the 
10 23 -rd zero. 



2A = 8 


2A = 10 


2A = 12 


0.50 


1.84 


4.51 


0.92 


3.32 


7.99 


1.24 


4.35 


10.2 


1.54 


5.35 


12.3 


1.77 


6.04 


13.7 



In the case of negative moments, the conjecture J_i(T) ~ A(logT/(2?r)) 1 
as T — s- co, due to Gonck [Go2 , suggests the negative second moment should be 
« 0.01808 near zero number 10 16 , w 0.01436 near zero number 10 20 , and w 0.01238 
near zero number 10 23 . These predictions are in good agreement with the values 
listed in Table El 

For 2A < —3, the behavior is much less predictable because, empirically, their 
sizes are determined by a few zeros where |£'(l/2 + ij n )\ is small. In fact, the 
particularly large fluctuations in the size of the negative sixth moment (2A = —6), 
near the 10 23 -rd zero in Table [SJ are essentially due to 8 zeros (out of 1.5 x 10 10 ) 
where |C'(l/2 + vy n )\ is equal to 0.002439, 0.002453, 0.002719, 0.002737, 0.003094, 
0.003108, 0.004365, and 0.004388. 

TABLE 6. Ratio j2.3| calculated with \B\ S3 10 9 , except for the first set, 
which uses the first 10 s zeros. The column "Zero" lists the approximate zero 
number near which block B is located. 



Zero 


2A = -2 


2A = -3 


2A = -4 


2A = -6 


10 s 


0.041129 





059025 


1.04212 


2935.6 


10 16 


0.018057 





030660 


0.55588 


1488.1 


10 20 


0.014341 





028403 


0.73586 


2873.2 


10 23 


0.012347 





022040 


0.41441 


1106.5 




0.012365 





022605 


0.43869 


1314.6 




0.012462 





037677 


2.76255 


63336 




0.012321 





021618 


0.42275 


1431.0 




0.012776 





178047 


59.6610 


9288238 




0.012326 





021062 


0.33853 


665.29 




0.012515 





052929 


7.46570 


412318 




0.012334 





022429 


0.56305 


4157.4 




0.012376 





025800 


0.81652 


5414.6 




0.012541 





089163 


21.5695 


2174342 




0.012411 





039415 


4.32860 


185114 




0.012329 





022729 


0.55154 


2723.6 




0.012386 





027487 


1.08706 


11563 




0.012605 





117993 


35.4067 


4686740 




0.012334 





021217 


0.33424 


538.73 



Starting with the investigations of |Od2] , several long-range correlations have 
been found experimentally in zeta function statistics. Such correlations are not 
present in random matrices, but do appear in some dynamical systems that for 
certain ranges are modeled by random matrices. So far all the zeta function corre- 
lations of this nature have been explained (at least numerically and heuristically) 
by relating them to known properties of the zeta function, such as explicit formulas 
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that relate primes to zeros. A natural question is whether such correlations arise 
among values of <^'(l/2 + ij n ). 

In order to detect correlations among values of |£'(l/2 + ij n )\, consider 

(2.4) S 2 (T,H,m):= ]T \C (1/2 + i ln )C (1/2 + i~f n+m )\ 4 . 

T<f n <T+H 

We computed this shifted moment function for various choices of m, T, and H. 
(We also considered similar sums with exponents other than 4, but for simplicity 
do not discuss them here.) Figure [3] presents some of our results near the 10 16 -th 
and 10 23 -rd zeros, and with H spanning about 10 7 zeros in both cases. The figure 
shows that correlations do exist and persist over long ranges. Also, the shape of 
S%(T, H,m) near the 10 16 -th zero is similar to that near the 10 23 -rd zero, except 
the former has higher peaks, and covers the range 3 < m < 222, as opposed to 
3 < m < 325, which suggests oscillations scale as 1/ log(T/27r). 

We remark the plot of S2 (T, H, m) in Figure [3] (right plot) is similar to a plot in 
[HOj of the shifted fourth moment of the zeta function on the critical line: 

rT+H 

(2.5) M(T,H;a):= \((l/2 + it)\ 2 \((l/2 + it + ia)\ 2 dt , 

Jt 

which we reproduce here in Figure [4] for the convenience of the reader. 



FIGURE 3 . Plots of S 2 (T, H, m)/S 2 (T,H,0) using 10 7 zeros near the 10 16 -th 
(left plot) and the 10 23 -rd zero (right plot) 




3 50 100 150 300 3 50 100 150 300 350 300 



To explain observed correlations, we numerically calculated the function: 



(2.6) 



f(T,H,x) 



53 C'(l/2 + *7n)e 

T<7„<T+i? 



27rinx 



which is related to long-range periodicities in £'(1/2 + *7n)- Assuming the RH, 
Fujii [Fu] supplied the following asymptotic formula in the case x = 0: 
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FIGURE 4. Plot of M(T,H,a)/M(T,H,0), with H rj 6.5 X 10 5 , near the 
10 23 -rd zero, drawn for a a multiple of 0.5. The dashed line is a sine kernel. 




(2.7) 

E C'(l/2 + ^)^log 2 £ + ( C o-l)£log£-( Cl+C o)£ +0 (TV 2 lo g ^ T ) , 

0<7„<T 



where Co = 0.5772 . . . (the Euler consta nt) and c\ = —0.0728 .... Empirical values 
of f(T,H,0) agree well with formula (2.7). For example, with H spanning 10 6 
zeros, we obtain f(T,H,0) = 21766088 - 14579i near the 10 20 -zero, and we obtain 
f(T, H, 0) = 25137126 + 61663i near the 10 23 -rd zero. But as x increases, f(T, H, x) 
experiences sharp spikes for certain x, as shown in Figure [5] which depicts the 
segment < x < 0.05 (in the remaining portion 0.05 < x < 1, the spikes get 
progressively denser). 

The sharp spikes in Figure[5]show the existence of long-range periodicities among 
values of ('(1/2 + ij n ). These spikes, as well as the correlations described above, 
are not unexpected. They can be demonstrated to follow from the properties of the 
zeta function, by estimating proper contour integrals. Such methods were used for 
continuous averages by Ingham Ingh] and even others before him, and for discrete 
averages over zeros by Gonek |Golj and Fujii |Fu[ IFu2j . The main step involves 
integration of C'(s) 2 /C(s), and estimates of such integrals. 

Applying such methods to (' (s) 2 e xsl ° s ^ /((s) suggests that the function 



f(T,H,x) 



E 

T<-f n <T+H 



C'(l/2 + i 7n )e 



In . T 

2^ l0g 2^ 



experiences large spikes at approximately x = log(fc)/ log(T/(27r)). For by a heuris- 
tic argument involving the (very) regular spacing of zeros one expects that j n in 
the definition of f(T, H, x) can be replaced by n without too much error (see |Od2j 
for a similar argument in the context of long-range correlations in zero spacings). 
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Therefore, f(T,H,x) should behave similarly to f(T,H,x)^\ In particular, we ex- 
pect the k-th spike in Figure [5]to occur at approximately log(fc)/ log(T/(27r)), and 
that agrees well with the evidence of the graphs. 

FIGURE 5. Plots of /(T,ff,x), defined in Ej& , using 10 6 zeros near the 
10 16 -th zero (upper left), 10 20 -rd zero (upper right), and 10 23 -rd zero (lower 
left). The lower right plot is another plot near the 10 23 -rd zero, except it uses 
a different set of 2 x 10 6 zeros. 



f(0) = 17372617 + 4917 i 



f(0| = 21766086 - 14579 i 



f(0) = 25137126 + 61663 i 



[(0| = 50220426 + 25141 



0.02 0.03 



0.05 0.00 



0.02 0.03 0.04 0.05 



3. Numerical methods 
As usual, define the rotated zeta function on the critical line by 

1/2 

(3.1) Z(t) = e»C(l/2 + it) , = ( ^4!^) ) ' 

The rotation factor e 10 ^ is chosen so that Z(t) is real. In our numerical experi- 
ments, t < 1.31 x 10 22 . 

Since |Z'(7„)| = |('(l/2 + «7„)|, it suffices to compute Z'(j n ). To do so, we used 
the numerical differentiation formula (Taylor expansion) 

(3.2) *(.)- '<' + »>-'<'-*> + «(«,», 

^Indeed, the plots in Figure [5] are almost unchanged if instead of plotting f(T,H,x) we plot 
f(T,H,x). 
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where the remainder term in (3.2 1 satisfies 

(3.3) \R(t,h)\< max l^-iM^ 
y ' t-h<H<t+h 6 

We chose h = 10 -5 , and approximated the derivative by 

(3.4) Z'it) « + -*(*-*) . 

To evaluate i/ (t) at individual points, we used a version of the Odlyzko-Schonhage 
algorithm OS implemented by the second author [Odlj . If the point- wise evalu- 
ations of Z(i + h) and Z(t — h) via this implementation are accurate to within 
±e each, then the approximation (3.4) is accurate to within ±(10 5 e + \R(t,h)\). 
Numerical tests suggested e is normally distributed with mean zero and standard 
deviation 10 -9 . Therefore, e is typically around 10 -9 . Also, varying the choice of 
h in (3.4 1 suggested the approximation is accurate to about 4 decimal digits with 
h = 10^ and t w 10 22 . 

In principle, our computations of £'(1/2 + i"f n ) can be made completely rigorous 
by carrying them out in sufficient precision. If one plans on calculating £'(l/2+i'y n ) 
with very high precision, however, it will likely be better to first derive a Riemann- 
Siegel type formula for Z'(t) itself, with explicit estimates for the remainder. Such 
a formula will be useful on its own as it can be be used to check other conjectures 
about C (1/2 + it). 



4. Conclusions 

Numerical data from high zeros of the zeta function generally agrees well with 
the asymptotic results that have been proved, as well as with several conjectures. 
There are some systematic differences between observed and expected distributions, 
but the discrepancies decline with growing heights. 

The results of this paper provide additional evidence for the speed of convergence 
of the zeta function to its asymptotic limits. They also demonstrate the importance 
of outliers, and thus the need to collect extensive data in order to obtain valid 
statistical results. The long-range correlations that have been found among values 
of the derivative of the zeta function at zeros can be explained by known analytic 
techniques. 
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